      subroutine aec8g2(coorr,coefr,
     & prmt,estif,emass,edamp,eload,num)
c .... coorr ---- nodal coordinate value
c .... coefr ---- nodal coef value
      implicit real*8 (a-h,o-z)
      dimension estif(24,24),elump(24),emass(24),
     & eload(24)
      dimension prmt(*),coef(3),coefr(8,3),
     & eexx(24),eeyy(24),eezz(24),eeyz(24),
     & eexz(24),eexy(24),eplast(24),
     & coorr(3,8),coor(3)
      common /raec8g2/ru(8,32),rv(8,32),rw(8,32),
     & cu(8,4),cv(8,4),cw(8,4)
c .... store shape functions and their partial derivatives
c .... for all integral points
      common /vaec8g2/rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      common /daec8g2/ refc(3,8),gaus(8),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(3),kdord(3),kvord(24,3)
c .... nnode ---- the number of nodes
c .... nrefc ---- the number of numerical integral points
c .... ndisp ---- the number of unknown functions
c .... nrefc ---- the number of reference coordinates
c .... nvar ---- the number of unknown varibles var
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      dimension dev(3),dep(3),e(3,3),d(3,3)
      dimension ddf(6),ds(6),p(4)
      external amohrc
      common /gpstr/ gpstr(7,1000000),str(7)
      nstr = 7
      pe=prmt(1)
      pv=prmt(2)
      fx=prmt(3)
      fy=prmt(4)
      Fg=prmt(5)
      p(1)=prmt(6)
      p(2)=prmt(7)
      p(4)=prmt(8)
      qnn=prmt(9)
      rou=prmt(10)
      wrou=prmt(11)
      grndwt_lv=prmt(12)
      alpha=prmt(13)
      edjt=prmt(14)
      time=prmt(15)
      dt=prmt(16)
      imate=prmt(17)+0.5
      ielem=prmt(18)+0.5
      nelem=prmt(19)+0.5
      it=prmt(20)+0.5
      nmate=prmt(21)+0.5
      itime=prmt(22)+0.5
      ityp=prmt(23)+0.5
      kq=int(qnn+0.5)
      fact = pe/(1.+pv)/(1.-pv*2)
      shear = (0.5-pv)
      p(3)=pv
      shear=shear*fact
      if (num.eq.1) call aec8g2i
c .... initialize the basic data
      do 10 i=1,nvar
      emass(i)=0.0
      eload(i)=0.0
      do 10 j=1,nvar
      estif(i,j)=0.0
10    continue
      do 999 igaus=1,ngaus
      call aec8g2t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr,
     & rctr,crtr,det,coefr)
c .... coordinate transfer from reference to original system
c .... rctr ---- Jacobi's matrix
c .... crtr ---- inverse matrix of Jacobi's matrix
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1,igaus)
      ry=refc(2,igaus)
      rz=refc(3,igaus)
      call eaec8g2(refc(1,igaus),coef,coorr,coefr,coefd)
c .... compute coef functions and their partial derivatives
      iu=(igaus-1)*4+1
      iv=(igaus-1)*4+1
      iw=(igaus-1)*4+1
      if (num.gt.1) goto 2
c .... the following is the shape function caculation
      call aec8g21(refc(1,igaus),ru(1,iu),rctr,crtr)
      call aec8g22(refc(1,igaus),rv(1,iv),rctr,crtr)
      call aec8g23(refc(1,igaus),rw(1,iw),rctr,crtr)
2     continue
c .... the following is the shape function transformation
c .... from reference coordinates to original coordinates
      call shapn(nrefc,ncoor,8,ru(1,iu),cu,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rv(1,iv),cv,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rw(1,iw),cw,crtr,1,4,4)
c .... the coef function transformation
c .... from reference coordinates to original coordinates
      call shapc(nrefc,ncoor,3,coefd,coefc,crtr,2,9,9)
      un=coef(1)
      vn=coef(2)
      wn=coef(3)
      weigh=det*gaus(igaus)
      do 100 i=1,24
      eexx(i) = 0.0
      eeyy(i) = 0.0
      eezz(i) = 0.0
      eeyz(i) = 0.0
      eexz(i) = 0.0
      eexy(i) = 0.0
      eplast(i) = 0.0
100   continue
      d(1,1) = +(1.-pv)*fact
      d(1,2) = +pv*fact
      d(1,3) = +pv*fact
      d(2,1) = +pv*fact
      d(2,2) = +(1.-pv)*fact
      d(2,3) = +pv*fact
      d(3,1) = +pv*fact
      d(3,2) = +pv*fact
      d(3,3) = +(1.-pv)*fact
      vol = 1.0
      if(z .le. grndwt_lv) then
        rou = wrou
      endif
      fz = Fg*rou
      if ((time .gt. 1.5D0*dt)
     +     .and. (num.eq.1) .and. (igaus.eq.1)) then
        open(17,file='filestr',form='unformatted',status='old')
        read(17) ((gpstr(i,j),i=1,nstr),j=1,nelem*ngaus)
        close(17)
      endif
      if (time.gt.1.5*dt) then
        do i=1,nstr
          str(i)=gpstr(i,(num-1)*ngaus+igaus)
        enddo
      else
        do i=1,nstr
          str(i)=0.0d0
        enddo
      endif
      e(1,1)=coefc(1,1)
      e(1,2)=coefc(1,2)
      e(1,3)=coefc(1,3)
      e(2,1)=coefc(2,1)
      e(2,2)=coefc(2,2)
      e(2,3)=coefc(2,3)
      e(3,1)=coefc(3,1)
      e(3,2)=coefc(3,2)
      e(3,3)=coefc(3,3)
      dev(1) = +e(1,1)
      dev(2) = +e(2,2)
      dev(3) = +e(3,3)
      dep(1) = e(2,3)+e(3,2)
      dep(2) = e(1,3)+e(3,1)
      dep(3) = e(1,2)+e(2,1)
      ds(1) = +d(1,1)*dev(1)+d(1,2)*dev(2)+d(1,3)*dev(3)
      ds(2) = +d(2,1)*dev(1)+d(2,2)*dev(2)+d(2,3)*dev(3)
      ds(3) = +d(3,1)*dev(1)+d(3,2)*dev(2)+d(3,3)*dev(3)
      ds(4) = shear*dep(1)
      ds(5) = shear*dep(2)
      ds(6) = shear*dep(3)
      ialpha=1
      str(1) = str(1)+ds(1)
      str(2) = str(2)+ds(2)
      str(3) = str(3)+ds(3)
      str(4) = str(4)+ds(4)
      str(5) = str(5)+ds(5)
      str(6) = str(6)+ds(6)
      call getyield(p,str,d,shear,a,ddf,kq,amohrc,ialpha)
      prag=amohrc(p,str)
      dln=+ddf(1)*dev(1)+ddf(4)*dep(1)+ddf(2)*dev(2)
     * +ddf(5)*dep(2)+ddf(3)*dev(3)+ddf(6)*dep(3)
      dln=dln+prag
      if (dln.gt.1.0d-6) then
        b=1.0/a
        ds(1) = +b*ddf(1)*ddf(1)*dev(1)+b*ddf(1)*ddf(2)
     1          *dev(2)+b*ddf(1)*ddf(3)*dev(3)
        ds(2) = +b*ddf(2)*ddf(1)*dev(1)+b*ddf(2)*ddf(2)
     1          *dev(2)+b*ddf(2)*ddf(3)*dev(3)
        ds(3) = +b*ddf(3)*ddf(1)*dev(1)+b*ddf(3)*ddf(2)
     1          *dev(2)+b*ddf(3)*ddf(3)*dev(3)
        ds(4) = +b*ddf(4)*ddf(4)*dep(1)
        ds(5) = +b*ddf(5)*ddf(5)*dep(2)
        ds(6) = +b*ddf(6)*ddf(6)*dep(3)
        str(1) = str(1)-ds(1)-ddf(1)*prag*b
        str(2) = str(2)-ds(2)-ddf(2)*prag*b
        str(3) = str(3)-ds(3)-ddf(3)*prag*b
        str(4) = str(4)-ds(4)-ddf(4)*prag*b
        str(5) = str(5)-ds(5)-ddf(5)*prag*b
        str(6) = str(6)-ds(6)-ddf(6)*prag*b
      else
        b=0.0d0
      end if
      ialpha=1
      call getyield(p,str,d,shear,a,ddf,kq,amohrc,ialpha)
      prag=amohrc(p,str)
      do 101 i=1,8
      iv=kvord(i,1)
      stif=+cu(i,2) 
      eexx(iv)=eexx(iv)+stif
101   continue
      do 102 i=1,8
      iv=kvord(i,2)
      stif=+cv(i,3) 
      eeyy(iv)=eeyy(iv)+stif
102   continue
      do 103 i=1,8
      iv=kvord(i,3)
      stif=+cw(i,4) 
      eezz(iv)=eezz(iv)+stif
103   continue
      do 104 i=1,8
      iv=kvord(i,2)
      stif=+cv(i,4) 
      eeyz(iv)=eeyz(iv)+stif
104   continue
      do 105 i=1,8
      iv=kvord(i,3)
      stif=+cw(i,3) 
      eeyz(iv)=eeyz(iv)+stif
105   continue
      do 106 i=1,8
      iv=kvord(i,1)
      stif=+cu(i,4) 
      eexz(iv)=eexz(iv)+stif
106   continue
      do 107 i=1,8
      iv=kvord(i,3)
      stif=+cw(i,2) 
      eexz(iv)=eexz(iv)+stif
107   continue
      do 108 i=1,8
      iv=kvord(i,1)
      stif=+cu(i,3) 
      eexy(iv)=eexy(iv)+stif
108   continue
      do 109 i=1,8
      iv=kvord(i,2)
      stif=+cv(i,2) 
      eexy(iv)=eexy(iv)+stif
109   continue
      do 110 i=1,8
      iv=kvord(i,1)
      stif=+cu(i,2)*ddf(1)
     & +cu(i,4)*ddf(5)
     & +cu(i,3)*ddf(6)
      eplast(iv)=eplast(iv)+stif
110   continue
      do 111 i=1,8
      iv=kvord(i,2)
      stif=+cv(i,3)*ddf(2)
     & +cv(i,4)*ddf(4)
     & +cv(i,2)*ddf(6)
      eplast(iv)=eplast(iv)+stif
111   continue
      do 112 i=1,8
      iv=kvord(i,3)
      stif=+cw(i,4)*ddf(3)
     & +cw(i,3)*ddf(4)
     & +cw(i,2)*ddf(5)
      eplast(iv)=eplast(iv)+stif
112   continue
c .... the following is the stiffness computation
      do 202 iv=1,24
      do 201 jv=1,24
      stif=+eexx(iv)*eexx(jv)*d(1,1)*vol
     & +eexx(iv)*eeyy(jv)*d(1,2)*vol
     & +eexx(iv)*eezz(jv)*d(1,3)*vol
     & +eeyy(iv)*eexx(jv)*d(2,1)*vol
     & +eeyy(iv)*eeyy(jv)*d(2,2)*vol
     & +eeyy(iv)*eezz(jv)*d(2,3)*vol
     & +eezz(iv)*eexx(jv)*d(3,1)*vol
     & +eezz(iv)*eeyy(jv)*d(3,2)*vol
     & +eezz(iv)*eezz(jv)*d(3,3)*vol
     & +eeyz(iv)*eeyz(jv)*shear*vol
     & +eexz(iv)*eexz(jv)*shear*vol
     & +eexy(iv)*eexy(jv)*shear*vol
     & -eplast(iv)*eplast(jv)*b*vol
      estif(iv,jv)=estif(iv,jv)+stif*weigh
201    continue
202    continue
c .... the following is the mass matrix computation
      stif=rou*vol
      elump(1)=stif*weigh
      stif=rou*vol
      elump(4)=stif*weigh
      stif=rou*vol
      elump(7)=stif*weigh
      stif=rou*vol
      elump(10)=stif*weigh
      stif=rou*vol
      elump(13)=stif*weigh
      stif=rou*vol
      elump(16)=stif*weigh
      stif=rou*vol
      elump(19)=stif*weigh
      stif=rou*vol
      elump(22)=stif*weigh
      stif=rou*vol
      elump(2)=stif*weigh
      stif=rou*vol
      elump(5)=stif*weigh
      stif=rou*vol
      elump(8)=stif*weigh
      stif=rou*vol
      elump(11)=stif*weigh
      stif=rou*vol
      elump(14)=stif*weigh
      stif=rou*vol
      elump(17)=stif*weigh
      stif=rou*vol
      elump(20)=stif*weigh
      stif=rou*vol
      elump(23)=stif*weigh
      stif=rou*vol
      elump(3)=stif*weigh
      stif=rou*vol
      elump(6)=stif*weigh
      stif=rou*vol
      elump(9)=stif*weigh
      stif=rou*vol
      elump(12)=stif*weigh
      stif=rou*vol
      elump(15)=stif*weigh
      stif=rou*vol
      elump(18)=stif*weigh
      stif=rou*vol
      elump(21)=stif*weigh
      stif=rou*vol
      elump(24)=stif*weigh
      do 301 i=1,nvard(1)
      iv = kvord(i,1)
      emass(iv)=emass(iv)+elump(iv)*cu(i,1)
301   continue
      do 302 i=1,nvard(2)
      iv = kvord(i,2)
      emass(iv)=emass(iv)+elump(iv)*cv(i,1)
302   continue
      do 303 i=1,nvard(3)
      iv = kvord(i,3)
      emass(iv)=emass(iv)+elump(iv)*cw(i,1)
303   continue
c .... the following is the load vector computation
      do 501 i=1,8
      iv=kvord(i,1)
      stif=+cu(i,1)*fx*vol
      eload(iv)=eload(iv)+stif*weigh
501   continue
      do 502 i=1,8
      iv=kvord(i,2)
      stif=+cv(i,1)*fy*vol
      eload(iv)=eload(iv)+stif*weigh
502   continue
      do 503 i=1,8
      iv=kvord(i,3)
      stif=+cw(i,1)*fz*vol
      eload(iv)=eload(iv)+stif*weigh
503   continue
      do 504 iv=1,24
      stif=-eexx(iv)*str(1)*vol
     & -eeyy(iv)*str(2)*vol
     & -eezz(iv)*str(3)*vol
     & -eeyz(iv)*str(4)*vol
     & -eexz(iv)*str(5)*vol
     & -eexy(iv)*str(6)*vol
     & +eplast(iv)*b*prag*vol
      eload(iv)=eload(iv)+stif*weigh
504   continue
999   continue
998   continue
      return
      end

      subroutine aec8g2i
      implicit real*8 (a-h,o-z)
      common /daec8g2/ refc(3,8),gaus(8),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(3),kdord(3),kvord(24,3)
c .... initial data
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      ngaus=  8
      ndisp=  3
      nrefc=  3
      ncoor=  3
      nvar = 24
      nnode=  8
      kdord(1)=1
      nvard(1)=8
      kvord(1,1)=1
      kvord(2,1)=4
      kvord(3,1)=7
      kvord(4,1)=10
      kvord(5,1)=13
      kvord(6,1)=16
      kvord(7,1)=19
      kvord(8,1)=22
      kdord(2)=1
      nvard(2)=8
      kvord(1,2)=2
      kvord(2,2)=5
      kvord(3,2)=8
      kvord(4,2)=11
      kvord(5,2)=14
      kvord(6,2)=17
      kvord(7,2)=20
      kvord(8,2)=23
      kdord(3)=1
      nvard(3)=8
      kvord(1,3)=3
      kvord(2,3)=6
      kvord(3,3)=9
      kvord(4,3)=12
      kvord(5,3)=15
      kvord(6,3)=18
      kvord(7,3)=21
      kvord(8,3)=24
      refc(1,1)=5.773502692e-001
      refc(2,1)=5.773502692e-001
      refc(3,1)=5.773502692e-001
      gaus(1)=1.000000000e+000
      refc(1,2)=5.773502692e-001
      refc(2,2)=5.773502692e-001
      refc(3,2)=-5.773502692e-001
      gaus(2)=1.000000000e+000
      refc(1,3)=5.773502692e-001
      refc(2,3)=-5.773502692e-001
      refc(3,3)=5.773502692e-001
      gaus(3)=1.000000000e+000
      refc(1,4)=5.773502692e-001
      refc(2,4)=-5.773502692e-001
      refc(3,4)=-5.773502692e-001
      gaus(4)=1.000000000e+000
      refc(1,5)=-5.773502692e-001
      refc(2,5)=5.773502692e-001
      refc(3,5)=5.773502692e-001
      gaus(5)=1.000000000e+000
      refc(1,6)=-5.773502692e-001
      refc(2,6)=5.773502692e-001
      refc(3,6)=-5.773502692e-001
      gaus(6)=1.000000000e+000
      refc(1,7)=-5.773502692e-001
      refc(2,7)=-5.773502692e-001
      refc(3,7)=5.773502692e-001
      gaus(7)=1.000000000e+000
      refc(1,8)=-5.773502692e-001
      refc(2,8)=-5.773502692e-001
      refc(3,8)=-5.773502692e-001
      gaus(8)=1.000000000e+000
      end


      subroutine aec8g2t(nnode,nrefc,ncoor,refc,coor,coorr,
     & rc,cr,det,coefr)
      implicit real*8 (a-h,o-z)
      dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10),
     & coorr(ncoor,nnode),coor(ncoor),coefr(nnode,*)
      call taec8g2(refc,coor,coorr,coefr,rc)
      n=nrefc
      m=n*2
      det = 1.0
      do 10 i=1,n
      do 10 j=1,n
      if (i.le.ncoor) a(i,j) = rc(i,j)
      if (i.gt.ncoor) a(i,j)=1.0
      a(i,n+j)=0.0
      if (i.eq.j) a(i,n+i) = 1.0
10    continue
c     write(*,*) 'a ='
c     do 21 i=1,n
c21   write(*,8) (a(i,j),j=1,m)
      do 400 i=1,n
      amax = 0.0
      l = 0
      do 50 j=i,n
      c = a(j,i)
      if (c.lt.0.0) c = -c
      if (c.le.amax) goto 50
      amax = c
      l = j
50    continue
      do 60 k=1,m
      c = a(l,k)
      a(l,k) = a(i,k)
      a(i,k) = c
60    continue
      c = a(i,i)
      det = c*det
      do 100 k=i+1,m
100   a(i,k) = a(i,k)/c
      do 300 j=1,n
      if (i.eq.j) goto 300
      do 200 k=i+1,m
200   a(j,k) = a(j,k)-a(i,k)*a(j,i)
c     write(*,*) 'i =',i,'  j =',j,'  a ='
c     do 11 ii=1,n
c11   write(*,8) (a(ii,jj),jj=1,m)
300   continue
400   continue
      do 500 i=1,nrefc
      do 500 j=1,ncoor
500   cr(i,j) = a(i,n+j)
c     write(*,*) 'a ='
c     do 22 i=1,n
c22   write(*,8) (a(i,j),j=1,m)
c     write(*,*) 'rc ='
c     do 24 i=1,ncoor
c24   write(*,8) (rc(i,j),j=1,nrefc)
c     write(*,*) 'cr ='
c     do 23 i=1,nrefc
c23   write(*,8) (cr(i,j),j=1,ncoor)
c     write(*,*) 'det =',det
      if (det.lt.0.0) det=-det
c     write(*,*) 'det =',det
8     format(1x,6f12.3)
      end

      subroutine aec8g21(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external faec8g21
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(faec8g21,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function faec8g21(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccaec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vaec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     faec8g21=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     faec8g21=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     faec8g21=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     faec8g21=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     faec8g21=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     faec8g21=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     faec8g21=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     faec8g21=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine aec8g22(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external faec8g22
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(faec8g22,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function faec8g22(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccaec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vaec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     faec8g22=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     faec8g22=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     faec8g22=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     faec8g22=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     faec8g22=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     faec8g22=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     faec8g22=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     faec8g22=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine aec8g23(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external faec8g23
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(faec8g23,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function faec8g23(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccaec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vaec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     faec8g23=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     faec8g23=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     faec8g23=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     faec8g23=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     faec8g23=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     faec8g23=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     faec8g23=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     faec8g23=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine taec8g2(refc,coor,coorr,coefr,rc)
c .... compute coordinate value and Jacobi's matrix rc
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(3),coor(3),coorr(3,8),coefr(8,3),rc(3,3)
      common /ccaec8g2/ x(8),y(8),z(8),un(8),vn(8),wn(8)
      external ftaec8g2
      do 100 n=1,8
      x(n)=coorr(1,n)
      y(n)=coorr(2,n)
      z(n)=coorr(3,n)
100   continue
      do 200 n=1,8
      un(n)=coefr(n,1)
      vn(n)=coefr(n,2)
      wn(n)=coefr(n,3)
200   continue
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dcoor(ftaec8g2,refc,coor,rc,3,3,1)
c .... coordinate value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function ftaec8g2(refc,n)
c .... coordinate transfer function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(3)
      common /ccaec8g2/ x(8),y(8),z(8),un(8),vn(8),wn(8)
      common /vaec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3) n
1     ftaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*x(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*x(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*x(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*x(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*x(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*x(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*x(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*x(8)
      goto 1000
2     ftaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*y(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*y(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*y(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*y(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*y(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*y(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*y(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*y(8)
      goto 1000
3     ftaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*z(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*z(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*z(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*z(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*z(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*z(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*z(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*z(8)
      goto 1000
1000  return
      end

      subroutine eaec8g2(refc,coef,coorr,coefr,coefd)
c .... compute coef value and their partial derivatives
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(3),coef(3),coorr(3,8),coefr(8,3),coefd(3,3)
      external feaec8g2
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dcoef(feaec8g2,refc,coef,coefd,3,3,2)
c .... coef value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function feaec8g2(refc,n)
c .... coef function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(3)
      common /ccaec8g2/ xa(8),ya(8),za(8),un(8),vn(8),wn(8)
      common /vaec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3) n
1     feaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*un(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*un(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*un(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*un(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*un(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*un(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*un(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*un(8)
      goto 1000
2     feaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*vn(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*vn(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*vn(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*vn(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*vn(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*vn(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*vn(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*vn(8)
      goto 1000
3     feaec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*wn(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*wn(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*wn(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*wn(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*wn(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*wn(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*wn(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*wn(8)
      goto 1000
1000  return
      end

